A STOCHASTIC ALGORITHM FOR PARAMETRIC SENSITIVITY 
IN SMOLUCHOWSKI'S COAGULATION EQUATION 
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Abstract. In this article a stochastic particle system approximation to the parametric sensi- 
tivity in the Smoluchowski coagulation equation is introduced. The parametric sensitivity is the 
derivative of the solution to the equation with respect to some parameter, where the coagulation 
kernel depends on this parameter. It is proved that the particle system converges weakly to the 
sensitivity as the number of particles N increases. A Monte Carlo algorithm is developed and vari- 
ance reduction techniques are applied. Numerical experiments are conducted for two kernels: the 
additive kernel and one which has been used for studying soot formation in a free molecular regime. 
It is shown empirically that the techniques for variance reduction are indeed very effective and that 
the order of convergence is 0(1/N). The algorithm is then compared to an algorithm based on a 
finite difference approximation to the sensitivity and it is found that the variance of the sensitivity 
estimators are considerably lower than that for the finite difference approach. Furthermore, two 
methods of establishing 'efficiency' are considered and the new algorithm is found to be significantly 
more efficient. 
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1. Introduction. Smoluchowski 's description of a coagulation process is made 
in terms of densities ^t(x) of particles of mass x = 1, 2, 3, . . . and takes the form of an 
infinite dimensional differential equation 

d 1 x — ^ °° 

^*( x ) = 2 XI K (yi x -y)L l t(y)vt{x-y) - Ht(x)^2K(x,y) fitiy). (1.1) 

The symmetric kernel K(x, y) appearing in this equation should be understood as 
giving the rate at which two particles of mass x and y coagulate. One gets an equiv- 
alent and more symmetric equation considering /it(-) as a measure on the set of 
non-negative integers and looking at the time evolution of observables of the form 
(/, /i t ) := ^2 X f(x)fj, t (x); moments are examples of such observables. In these terms, 
equation Ijl.ip takes the form 



(/,!"*) = (/, Mo) + \J a f £ {f{x + y)-f{x)-f{y)}K(x,y)tx s {x) i ji s {y)\ ds. 

(1.2) 

The basic problem we address is to derive a numerical scheme to understand how 
the solution to this equation depends on possible parameters in the kernel. We shall 
write K\ to indicate that K depends on some d-dimensional parameter A, and shall 
write for the solution of equation (|1.2p . Formally differentiating this equation with 
respect to A and setting er A = dx^t we g et 

(,/> t A ) = (,/> A ) + - / {f(x + y)-f(x)-f{y)}K x {x,y)ri{x)a*(y))ds 

+ / (E {f( x + y)-f( x )-f(y)} K x(x>y)^(x)^(y)) a*. 

(1.3) 

K' x is here the derivative of K\ with resrJect to A. Section [2] presents an algorithm 
which simulates the sensitivity ct a very accurately and in an efficient way. 
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There are two main motivations for performing sensitivity analysis. The first 
is for solving inverse problems. If some particle system is governed by a partial 
differential equation which in turn is dependent upon some unknown parameter, it 
is desirable to find this parameter. This can be achieved by choosing the parameter 
value which minimises some residual which is a function of experimentally realised 
quantities and its computational analogue, which varies with the parameter. The 
minimisation procedure often uses a gradient search, thus the value of computing 
parametric derivatives is apparent. Secondly, in considering a scientific model, we 
often wish to consider the smallest model which reasonably fits the data, in which 
case sensitivity analysis can be performed to discard parameters with small sensitivity. 

Whilst the usual tools of solving differential equations (and their associated nu- 
merical schemes) are badly adapted to the above infinite dimensional framework, the 
stochastic approach of interaction particle systems (basically Markov chains) can be 
used efficiently, in this setting, as Marcus in pQ, and later Lushnikov in [2], first re- 
alised. We follow their approach and give a stochastic particle approximation of the 
sensitivity af. 

Before running any simulation, one should investigate the well-posedness of equa- 
tion l|1.3|) : if it had more than one solution it would be unclear what solution a nu- 
merical scheme approximates. The most general answers to this theoretical question 
for Smoluchowski equation were given by Jeon in [3] and Norris in [4], under a growth 
assumption on the interaction kernel and a moment condition on the initial condition 
fio. Surprisingly enough, the existence and uniqueness problem for the sensitivity was 
only solved recently, by Bailleul [5], using methods developed by Kolokoltsov [6]. The 
algorithm developed in this article is the numerical counterpart of this theoretical 
wor!0. 

Three approaches to the simulation of the sensitivity by systems of particles 
have mainly been used up to now. The first uses weighted particles, as illustrated 
by Vikhansky and Kraft [7]. They approximate the family of solutions {^t} \ by 
Marcus-Lushnikov processes J2n^o Wn (^ > ^)^„(t) where the dependence on A is en- 
tirely put on the weights w n (t; A). A heuristic argument imposes to their derivative 
to satisfy a kind of Markov evolution rule. Despite its (numerically verified) conver- 
gence this approach essentially has the same speed of convergence and variance as the 
Marcus-Lushnikov process. Further, the paper does not any information regarding 
computation run times. 

The second approach considers adjoint sensitivity [8]. A backward partial dif- 
ferential equation is used rather than a forward one, as used in most other methods. 
The advantage of this method is that sensitivity for any parameter value is immediate 
once the computation have been done whereas using the forward equation requires 
explicit calculation for each parameter value. The disadvantage is that one can only 
calculate the sensitivities for a particular functional of the particle ensemble. 

In the third approach, devised by the authors with J. R. Norris in the forthcoming 
article [9], the sensitivity a x is approximated by the ratio (fi x+dX ' N — fi X,N )/SX, 
where ^ X+SX ' N and ^ X ' N are two Markus-Lushnikov processes corresponding to close 
parameters, coupled so as to minimise the difference of their random fluctuations 
around [i x+SX and fi x . This approach leads to a massive decrease of variance but does 
not improve the speed of convergence of the algorithm. 



1 Consult this article for conditions under which existence and uniqueness of a solution to the 
sensitivity equation l|1.3|l holds. 
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The algorithm we propose improves the variance of the sensitivity estimator and 
requires a much smaller number of particles to converge. This is described in section 
[21 The reader who is not interested in mathematical details can skip sections 12.11 and 
12.21 where it is proven that the particle system introduced in section [2] converges to 
the sensitivity. Section [3] presents the algorithm we have used to obtain the numerical 
results of section 01 

Notation. We shall prove convergence of the particle system in a general setting 
where masses of particles can take any positive real value. The densities of particles 
will then be represented by non-negative measures fit and all sums will be replaced by 
integrals. In this framework we shall write (/, /i) for J f(x)fi(dx) and Smoluchowski's 
equation l|1.2p will be written 

(/, IH) = (/, Mo) + \ J j {f(x + y)- f{x) - f(y)} K(x, y) ^ s (dx) fi s (dy) ds. 
We shall formally write it as 

(4 = ±KM^t)- (1.4) 
In the same way, we shall write formally equation (|1.3p for the sensitivity as 

= K\ a^) + \k' x (Mt , fJ>t) ■ (1.5) 

The integral notation is adopted from now on. 

2. Markov chain approximation. It is probably fair to say that although the 
Smoluchowski equation l|1.2p is a deterministic evolution equation it should primar- 
ily be thought of as a deterministic large scale picture of a stochastic mesoscopic 
dynamics. Indeed, Smoluchowski obtained his equation from a representation of the 
coagulation process using 'particles' moving according to Brownian trajectories whose 
diffusivity depends on their mass and coagulate when they are close to each other. 
As explained in the article [10] of Chandrasekhar, section 6 of chapter III, in a region 
of space where the coagulating particles are well mixed, one can forget about their 
spatial location and obtain a mean-field evolution for their mass distribution. This 
mean-field picture is provided by Smoluchowski equation. Given in its simple form 
(|l.ip . it is not clear at first sight how one should simulate a solution to this infinite 
dimensional differential system. 

The approach developed by Marcus in his seminal paper [I] in a sense comes 
back to the primary stochastic description of the coagulation phenomenon and relies 
on the intuitive content of Smoluchowski equation. Two particles of masses x and 
y coagulate at rate K(x,y) to create a new particle of mass x + y: The particles x 
and y are removed from the system and the particle x + y added. This motivated 
Marcus, and later Lushnikov, to represent a particle of mass x by a Dirac mass S x at x 
and to introduce a strong Markov jump process on the space of discrete non-negative 

measures with the following simple dynamics. Denote by = — 8 Xi its initial 

i 

state and by ^ its state at time t. Associate to each pair 1 ^ i < j ^ N independent 

K(xi, Xj) 

exponential random times Ty with parameter — — — — - and set 



T = min{T„ ; 1 < i < j N}. 
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The process remains constant on the time interval [0, T), and if T = T pq it has a jump 
— (S Xp+Xq — S Xp — S Xq ) at time T. The dynamics then starts afresh. Note that the 

new measure at time T is still non-negative, and that the above description leads to 
a mean jump of the process during a time interval [t, t + St] equal td3 

£*£(W - S x - S xl ) K(x, y) tf(x) /if (a/) 

x,x' 

up to terms of order jj and o(St). This property makes it clear that the process 
converges to a solution of the Smoluchowski equation as TV goes to infinity (under 
proper conditions), a fact which was used for simulation purposes long before it was 
proved under general conditions in [4]. 

Following the heuristic approach of Marcus and Lushnikov, we are going to give 
in the next section a particle description of the sensitivity equation 

a t A = K x (rf, a f A ) + \k' x (^1 rf). (2.1) 

To that end, introduce the notation K' + := K' V and K'_ := K' A (dropping the 
index A for it will be fixed), and write, for a signed measure a, 

v = rjl -^ 1 >Q - Ml <0 =: ct + - (7~, 

d\a\ d | cr | 

Using this notation, re-write equation l|2.ip as 

*t - K = (^A(MtW) + \K'M^t)) - (#a(m*W) + \k'_{^,^)) (2-2) 

This equation will motivate the introduction of the Markov chain described in the 
next section. 

Notation. Given three non-negative measures fi, er + , a~ on (0, oo) we shall adopt the 
notation \x © a + @ <j~ to denote the Revalued measure on (0, oo) 3 . It will clarify the 
notation to denote by x © y © z the point of R 3 with co-ordinates x, y and z. Given 
non-negative functions /, g, h on (0, oo) set 

(/ © g © h, n © cr + © a') := (/, © (<?, a+) © (h, a~). 

As we shall simulate both fi t and (er^",cr t - ) at the same time, our approximating 
Markov chain will take values in the set 

TV := {/! © er + © <j~ ; fi, a + , a~ non-negative discrete measures on (0, oo)}. 

2.1. Chain, generator. In the same way as the right hand side of Smoluchowski 
equation (|1.4p can be interpreted as the coagulation of particles of fj, t of mass x and y 
at rate K(x,y), we are going to follow what equation (12.21) suggest and interpret the 
term K({i t , &t) appearing there as the coagulation of a particle in /it of mass x with a 
particle in of mass y at rate K(x, y). Note that this leads to a jump S x + y — 8 X — S y 
of cr + which could transform the non-negative measure af into a signed measure, as 
the term 5 X does not necessarily appear inside af (while S y does) . We shall take care 
of this by adding S x to the negative part <j^~ of at instead of subtracting it from erf ; 



2 fi^ denotes the state of the process at time t. 
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as we are only interested in the difference er t + — <r^(= a t ) this has no consequence. 
Note also that the particle 8 X from fi t used in that coagulation event will not be 
removed from /it. Similar interpretations of the terms K(fi t ,o'^) and ^K' ± (^i t , fi t ) 
lead us to define the following Markov chain 8t = X t ffi Y t ffi Z t on A/\ Denote by 

©o = ( ^2 ^xj) ® ( ^ ^j/fej ® f fize) its starting point. 

i=l..m k—l..p £—l..q 

2.1.1. Dynamics. Associate to each pair 

• 1 < i < j < m, exponential random variables Rij,Sij and Ty with respective 
parameters K(xi,Xj) and K' + (xi,Xj) and K'_{xi,Xj), 

• {i,k) G x an exponential random variable f/jfc with parameter 

K{x u y k ), 

• (i, G [l,m] x [l,g]| an exponential random variable Vu with parameter 

K(xi, zi). 

All these random variables are supposed to be independent. Denoting by W the first 
event happening in the system 

W = mm{R lJ ,S lJ ,T lJ ,U lk , Vu ; 1 < i < j < m, fe G [l,p], * G 

the jump AO of the Markov chain depends on which of these exponential clocks rings 
first. For future reference, the different types of events that can happen are numbered. 
If 



w = 


R-ij 3 


then 


A9 = 


w = 




then 


A6 = 


w = 


T- 


then 


A6 = 


w = 


Uik, 


then 


A6 = 


w = 




then 


A6 = 



) © © (event type 

ffi 8 Xi+Xj ffi (8 Xi + 5 Xj ) (event type 

ffi (8 Xi + S X j) © S Xi+Xj (event type 

ffi \S Xz +y k - S yk ) ffi 8 Xi (event type 

ffi S Xi ffi (S Xi+Zt - 8 ze ) (event type 



) 

1 + ) 

1") 
2+) 

2") 



The process 0( will be constant on the time interval [0, W) and have jump A6 at 
time W . The dynamics then starts afresh. 

Remark. It is clear from this description that for any function ip satisfying the 
relation ip(a + b) ^ <p(a) — tp(b) for any a, b > 0, the function (ip, Yt + Zt) increases 
with time. This fact is useful for the convergence result stated in theorem \2.2[ 

Given any positive integer N, define jjQt as the element jjXt ffi jjYt ffi jf%t of 
Af, and set 

c N N 

Note that the first component of 9^ is the usual Marcus-Lushnikov process. Set 
— af' N — &t' N ■ We are going to prove in theorem 12.21 that <r^ converges in law 
to the sensitivity at- Those who do not care about the mathematical details of such a 
statement can skip the remaining of this section and go to section^ 

2.1.2. Generator. The analytic description of the Markov chain {6^)4^0 in 
terms of its generator will be useful in proving theorem 12.21 Given a non-negative 
measure fx of the form define the rescaled counting measure on ordered pairs 

of masses of distinct particles by 

?{A x A') := »(A) fJ-(A') - ±n(A n A'), 
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and define the measure G^ N ^(fj,) and the operator P( n \(jl) setting for any measurable 
bounded function / 

(/,GW( M ))=i J{f(x + x')-f(x)-f(x')}K(x,x')Jl(dx,dx / ) 

(/,P m W)4 [{f(x + x')-f(x)-f(x')} 2 K(x,x')jl(dx,dx'). 

Given x > and a non-negative measure 7 on we shall write K(x,j) for the 
integral J K(x, y)j(dy). 

Denote by the generator of the process {®t'} 0<t<T ] f° r anv bounded mea- 

surable functions /, g, h on (0, 00) the K 3 -valued process 



Mi 



f,9,h;N 



{f®g®h,Q?) - (f®g®h,e^)-J Q (f®g® /i.H^ef)) ds 



is a martingale (with respect to its natural filtration) . For a measure /x of the form 
jf $xi and = fj, ® <t + © cr~ ejVwe have 



(/© 5 ©/i,H (w) (e) 

(/,G W W) © 

{\J{s( x + X ') K '+ + (g(x) +g(x'))K'_(x, x')}jl(dx, dx') 
+ J {(g(x + y) ~ g(y))K(x,y)a + (dy) + g(x)K(x,a~)}n(dx)^ 
i I \h{x + x')K'_{x,x') + (h(x) + h(x'))K' + (x,x')}ji(dx,dx') 
(h(x + z) - h(z))K(x, z)er~(dz) + h(x)K(x, a + )^{dx) 



(2.3) 

Compare this formula with the description of the dynamics given in the section 

(i) Event {W = R tj } corresponds to the term (/, G (Ar) (/x)) © © 0; 

(ii) Event {W = S t j} corresponds to the term \ J © g(x + y) © (h(x) + 
h(y)) K' + (x,y)Jl(dx,dy); a similar term corresponds to the event {W = Tij}; 

(hi) Event {W — Uik} corresponds to the term J{0 © (g{x + z) — g(z)) ® 
h(x)K(x, z)a + (dz)} fi(dx); a similar term corresponds to the event {W = Vu}. 
The sum of all these terms gives (/ © g ® h, H("'(9)). 

Following a classical approach, the study of martingales of the form M/' 9 ' /i ' N will 
be our main tool in the proof of the convergence theorem. The explicit expression of 
the bracket of M^ 9 ' h ; N will be useful in that task. We have 

(Mt** ; N ) t = 1 jf* (/ © g © h, Q <*> (6f )) ds, 
where Q Ar (0) is characterised on measures of the form {■^^28 Xi ) © c + © cr~ by 



A STOCHASTIC ALGORITHM FOR PARAMETRIC SENSITIVITY 



7 



the formula 

(/,P W (M)) © 
\ f {g(x + x' fK' + {x, x') + (g{x) + g{x')f K'_(x, x')]jl(dx, dx') 

[g(x + y)- g{y)) 2 K(x, y) a+(dy) + g{xf K(x, cr")} fj,(dx) 

i j \h(x + x') 2 Ki{x,x') + (h{x) + h{x')) 2 K' + {x,x')}'il{dx,dx') 

(h(x + z) - h(z)) 2 K(x, z)a~(dz) + h{xf K(x, a + )}n(dx) 



2.2. Convergence theorem. Denote by U a bounded open set of some M. d 
indexing the family K\ of kernels. Let ip : (0, oo) — » R+ be a sublinear function: 
ip(sx) ^ s<p(x) for any s > and x 6 (0,oc); such a function is also subadditive: 
ip(x + y) < (p(x) + ip(y), for any x, y € (0, oo). We shall suppose that the interaction 
kernels K \ satisfy the growth condition 

K x (x,y) sC (p(x)ip(y) 

for any x, y € (0, oo), AgW, and that the initial condition of Smoluchowski equation 
(|1.2|) (or better its 'continuous mass version') satisfies the moment condition 

V?(x) 4+e /! (da;) < oo (2.4) 



for some (small) e > 0. We shall suppose in theorem 12.21 that <p 2 is sub-additive; 
together with the above moment condition l|2.4p on fio this implies that Smoluchowski 
equation has a unique strong solutiorH, defined for all non-negative times. 

The following norm was used on the space M. \ of signed Borel measures /i such 
that ||/i||i := (</?, | /i |) < oo, in the article [5] where the following key result about 
sensitivity is proved. 

Theorem 2.1. Assume the moment condition l|2.4p and that K\(x,y) and 
\K' x (x,y)\ are both bounded above by tp(x)<p(y) for any x,y. Then the map (t,X) E 
[0, oo) x U i— > fit G (.Mi, ||.||i), is a C 1 function and its derivative satisfies the 
following equation for any bounded measurable function /fl). 



(/,<#) = + J J {f}{x,y)K x (x,y)Hs(dx)cr*(dy)ds 

{f}(x, y)K' x (x, y)ri(dx)^(dy)ds 



The function <r x is the only (Mi, ||.||i)-va/raed solution of this equation. 

We shall consider here a weaker topology than the || • ||i-topology. We shall equip 
the space R^ 3 with the ^-distance: \\x®y®z-x'®y'®z'\\ := \x-x'\ + \y-y'\ + \z-z'\. 



3 In the sense denned in [3]. 

4 We write here {f}(x,y) for f(x + y) - f(x) - f{y). 
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Write M® 3 for the set of non-negative R® 3 -valued measures on R+, and let d be any 
distance on M® 3 metrising weak convergence: {0„}„^o converges to 0oo iff for any 
bounded continuous functions /, g, h on R+, we have (/©<?ffi/i, 9„) — * (/©5©/i, ©oo) • 
The space (.M® 3 , d) is a Polish space with M as a dense subset. 

Fix a positive time T . We shall state our convergence theorem in the functional 
setting V([0,T],(M m ,d)) of cadlag paths from [0,T] to (M m ,d). This space will 
be equipped with its Skorokhod topology, for which we refer the reader to the books 
[TT] or [32] of Billingsley and Pollard. Last, we shall denote by d any distance on 
the set of all non-negative Borel measures on (0, oo) metrising the following notion 
of convergent: {fJ- n }n^o converges to iff we have (/, yu„) — > (/, Moo) for any 
bounded continuous measurable function / with bounded support 

The starting point 9^ of Q N will be of the form j^X^ © ±Y N © jjZfi for some 
non-negative integer-valued finite measures , , Z$ on (0,oo). To shorten the 
notation we shall denote by 

the process starting from Qq constructed in section 12.11 and corresponding to a given 
parameter A. 

We shall suppose that the function <p> controlling the kernels K\ satisfies identity 
(|2.5|) below. As noted in the remark on page this hypothesis implies that the 
function (ip, of' +cr^' ) increases with time; this fact will enable us to control 0^. 
Note that this hypothesis is weaker than requiring that ip be increasing. 

Theorem 2.2 (Convergence of the particle system). Let K\{-, ■) : W + x W + — > 
[0, +oo) be a family of symmetric kernels indexed by A G U. We suppose the map 
(A ; x, x') i— > K\(x, x') continuous and differ entiable with respect to X, with a derivative 
K' x (x, x') continuous with respect to (x, x'). Let ip 1 be a subadditive function whose 
square is also subadditive. Assume that 

tp(a + b) ^ fip) — <£>(&), for any positive a, 6, (2-5) 



VXeU,Vx,x',yeR* + , K x (x,x') ^ <p(x)tp(x'), 

i i (2-6) 

\K x (x,y)\ < tp(x)(p(y), 



K ^ x ') and K '^') _ o (2.7) 

ip(x) <fi(x') (p(x)ip(x') x+x'^oc 

Fix X G U and write <d N for the corresponding process inAf, started from ^ ®o~q' N © 
(Tq' . Suppose that /io satisfies the moment condition (|2.4p for some (small) e, that 

do(<pi<$,<piio)->0, (2-8) 

and that there exists a positive constant C bigger than {ip 2 , (Xq) and Up, Oq~' +a" ' N ) 
for any N ^ 1 . 



5 This notion of convergence, usually called vague convergence, is weaker than weak convergence. 
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Then the sequence of the laws of the processes is tight and any (random) weak 

0<t<T' 



limit is almost surely of the form {fit © © a t °°}o<t<T' m ^h 



+.00 — ,00 



Proof. The following estimate is essential in controlling the behaviour of the 
processes a + ' N and <r~' N . 

Lemma 2.3. There exists a positive constant C\ such that 



I +,N . -,N 

sup (if, a t +a t ' 

0<t<T 



First decompose Up, af' N + a t as the sum of a martingale {M t } Q<t<T and a 
finite variation term: 



{<p(x + x') + <p(x) + ip(x')} K'(x, x') Jif(dx, Ax') 
+ J { V (x + y)-rty) + rtx)}K(x,y)^(dx)(<jf< N + a-< N )(dy)^j d 
From j2j| we have for each N ^ 1 and t <E [0, T] 
(<p,a+- N + at' N ) <C + M t + [ f 2{ip(x) +i p(x')}cp(xMx')^(dx)^(dx')d 

U*t'» + a 7 > N )ds. 



s 



This upper bound is simplified using the subadditivity of <p and <p 2 from which we 
havcH 

(<p,tf)<{<p,IJ$)<C and C. 
This gives a Gronwall-type inequality 

(tp, ar+' N + a;- N ) ^C + M t + AC 2 T + 2C f Up, a+' N + a~' N ) ds 



whose mean version gives a constant C\ such that E Up, ut' N + cr t :N ) C± for any 
^ t ^ T. We get the statement of the lemma recalling that hypothesis 1|2.5|1 implies 
that the function 1 1— > Up, a^' + a^' ) is increasing. 
Given e > define the compact subset 

K t = |/i © cr + ® o~ e M m ; max{(<p,n), (^,cr + ), ((p,<x~)} < -} c X® 3 , 



3 Since f^lwe have /Xq ) ^ (ip 2 ,/^) ^ C. 
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and denote by V N the law of 9 W on X>([0, T], (M m , dfj . 

Corollary 2.4 (Compactness). Given r\ > 0, there exists e > such that 

¥ N (D([0,T],K e )' S j >l-r). 



Now let /, g, h be bounded measurable functions on (0, oo) no greater than 1. By 
lemma 12.31 we have for all s < t 



E 



(f®g®h,H.W(Q?))\\d8 
) + 2 



and 



2C 2 (t - i 
«S C 2 {t-s) 



E ^ M f,9,h;N^_( M f,g,h;N^ 



3C 2 
2 



— + 2C( V ,a r 



dr 



< —E 

^ N 



AC 2 

~w 



1 

N 



(f®g®h,Q {N) (Q?))\\d S 
[C 2 +AC 2 



2E 



AC Uprf 



C(<p, a r 



-,N\ 



(]>■ 



where C2 is a positive constant depending only on C, So, by Doob's L 2 -inequality, 
we have 



3 



sup \\{f®g®h,9? -ef)| 



(2.9) 



for some positive constant C3 depending only on C. It is then a standard fact that the 
equicontinuity inequality l|2.9p together with corollary on compactness enable the use 

of Jakubowski's criterioifl; so the sequence of laws of Q N inX>f[0,T], (M® 3 ,d)^j has a 

convergent subsequence. Denote by 8°° = 0cr + '°° 0O- - ' 00 any limit point. Taking 
a subsequence and changing the probability space if necessary we can suppose without 
loss of generality that 6 converges almost surely to 6°° in T>([0, T], (M m , d)Y As 

Q N makes jumps of size at most jjj, in the total variation distance, the limit process 
is a continuous process from [0, T] to (.M ffi3 ,d). 

It is proved in [4] that under conditions l|2.8|) and (|2 .4[) the process is almost 
surely equal to the unique strong solution /i. of Smoluchowski equation, and that we 
have almost surely sup do(<pti^ , f^ s ) — » 0, as N goes to 00. 

To prove that (T + '°° — cr. - ' 00 is equal to the unique solution of equation l|1.5p it 
suffices to prove that it satisfies this equation for any bounded measurable function g 



7 See for instance Dawson's lecture notes [13] 
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with compact support, as a straightforward limit argument will give it for any bounded 
measurable function. We shall suppose without loss of generality that a ' — aZ' = 
0. We shall adopt the notation 



oZ — er ' , 



-,N 



-,N 



and 



The conclusion of lemma 12.31 can now be re- written as E 



sup \}p,\(r t 

0<t<T 



N I 



It can be seen from expression (|2.3p for that the real-valued process 



B 



g;N 



: {g(x + x') - g(x') - g(x)} K'(x, x') ^ (dx, Ax') 



(2.10) 



+ j {g(x + v)- g(y) - g(x)} k(x, y) (dx) of (dy)J ds 

is a martingale with previsible increasing process 

(BO = 1 ^ (7 i{g{ x + x') g(x') g(x)} 2 K'(x, x') Jif (dx, dx') 

+ J {g(x + y) - g(y) - g(x)f K(x, y) ^ (dx) a?(dyf) ds 

Using lemma [231 together with the almost sure inequality (<£>, (i^) < C, it is seen that 
E( x B 9 ' ,N ) T converges to as N goes to oo. So, to show that a°° satisfies equation 
(|1.5|) . it is sufficient to prove that the two integrals inside the right hand side of 
equation l|2.10p converge almost surely to 



and 



^{g(x + x') - g(x') - g(x)} K'(x,x') ^ s (dx) n 8 (dx') 



J {g(x + y)- g(y) - g(x)} K(x, y) ^ s (dx) af(dy) 



(2.11) 



respectively, and that we have uniform bounds on them so that dominated convergence 
under the time integral can be used. The convergence of the first integral was proved 
in [4] using hypotheses (|2.6j) and (|2.7p . with K in place of K'\ the same argument 
applies here. This integral is bounded above by | ||fli|ooC 2 , uniformly in s £ [0, T] and 
N ^ 1. 

Given S £ (0, oo], the function <p s (x) = f(x)l x ^s is subadditive. It comes from 
Fatou's lemma that the inequality 



E 



sup (^Jcr^l 
0<t<T 



holds for any 5 £ (0, oo]. So, to any lu £ f2 one can associate a positive constant 
m(S ; u>) such that we have 

(/,|0)|) ^ {<p s ,\aF(u)\)^m{S;u>) 
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on the time interval [0, T]. One can choose this constant m{8 ; oj) so that it converges 
to as 5 decreases to 0. Taking w in a subset fl\ of f2 of probability 1, for which 
e N (cj) converges to e°°(w) in X>([0,T], (.M® 3 ,d)), we get that 

(<p S ,\o?(u,)\)^{<p 8 ,\*${u)\) 

is arbitrarily small provided 5 is small enough, and bounded above uniformly in t G 
[0,T],N^ 1 and S € (0,oo]. 

Proceed now as in [4] and write K as the sum of a kernel K\ with compact support 
and a kernel K 2 with support in 



UF 2 UF 3 :={{x,y);x^6}u{(x,y);y^S}U {(x,y) ; max{x,y} > i}. 



There is no problem in justifying the convergence of the integral in (|2.11j) correspond- 
ing to K x . For K 2 write, with {g}(x, y) := g(x + y) - g(x) - g(y), 



{g}(x, y) K 2 (x, y) ( ^ {Ax) af (dy) - p a (dx) a s °°(dy) 



N i 



{g}(x, y) K 2 (x, y) (/if - fif) (dx) af (dy) 
{g}(x, y) K 2 (x, y) (i s (dx) (<rf - af ) (dy) 



and deal with each term of the right hand side separately. The first term is bounded 
above by do(wf (a>), t^Ms) (f, |of (w)|), up to a multiplicative constant. As the first 
factor converges to (and is no greater than 2C) while the second is uniformly 
bounded above, one can apply dominated convergence in the corresponding integral 
with respect to s. To deal with the second term, use the pointwise bound^J 

\\K 2 l Flf i s ©af HI,, < j s C(<p, |af 
ll^l^^eafMH^C^^IafKo;)), 
||if 2 l F 3 Ms ©af(o/)|| < (/,/i s )(v>f 

where 75 = maxi ; (x,y) e F3J converges to as S decreases to 0. As 

Up, |crf (w)|) is uniformly bounded above by a constant, and both Up s , \af (w)|) and 
Up 8 , /i s ) can be made arbitrarily small for small enough S, we have enough control to 
apply dominated convergence. □ 

3. Algorithm. We describe in this section the algorithm used to simulate the 
particle system studied above; the numerical results are to be found in section |4j 
Two points of computational interest are first put forward in sections 13 . 1 1 and f3T2| the 
algorithm itself is described in section 13.31 



3.1. Coupling. The basic algorithm to simulate the sensitivity a t is given by 
the dynamics of the process Q N described in section 12.11 A fresh look at it reveals 
a potential computational drawback of this approach: It is seen from the explicit 
expression l|2.3p of the generator of 6^ that the mean number of particles inside 
<j N satisfies a Gronwall-type inequality, which implies an exponential growth of this 



|o denotes total variation norm. 
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quantity. One should see in this exponential growth of the number of particles a good 
feature for the approximation qualities of our estimator cr^ of a t , especially regarding 
accuracy and variance. This should be opposed to what happens for the weighed 
and coupled particles systems described in the introduction, for which the number of 
particles in the system decreases with tim^|. 

As an exponential growth of the quantity of information to consider is non- 
desirable for simulations, three kinds of tricks are used in order to reduce it. 

(i) Cancellation. As we are only interested in the difference — cr^' N any 
particle which appears in both particle systems will be removed from both of them. 

(ii) Coupling. A particle S x of fi N coagulates with any particle of a^' N at rate 
jtK (x, cr^~' N ) = -A- J K(x,y)a^' N (dy); it also coagulates with any particle of a^' N at 
rate jfK(x, &t' N ) ■ This particle is thus used in both systems at rate jjK(x, ot' N ) A 
K(x, <jf' ), in which case a cancellation removes the particles S x added to a t ' N and 
a +' N _ T/ ms operation leaves the total number of particles in <r N constant. The rest 
of the time 8 X is used in only one of the systems. 

(hi) Re-sampling. A more drastic control of the number of particles in a N can 
be obtained using re-sampling. Let M and m be two integers depending on N, with 
m < M . Each time ap N or a^' N has M particles, replace it by an iid sample of itself 
of size to; this way the total number of particles in a N remains no greater than 2M. 

3.2. Majorant kernel. In order to treat information in a computationally effi- 
cient way, we have organized the data using tree structures. The use of a majorant 
kernel with a simple algebraic structure together with an acceptance/rejection step 
lead to an efficient updating of the data tree. 

The choice of a majorant kernel K(-, •) is made so that K is symmetric, no less 
than K and has the form 

K(xi,Xj) = ^2K p {x l ,x j ) :=^2f (xi)g (xj) (3.1) 



for P in a finite set of indices [14] . This form of kernel leads to simple generation of 
probabilities of the form 

K(xj,Xj) = Eq/J^^lguM fpjxj) gp{x 3 ) 

where a and b run in possibly different finite sets of indices. Identity f|3.2|) corresponds 
to choosing first an index (3 according to the probability specified by the first term of 
the right hand side and then choosing each particle Xi,Xj separately. The choice of a 
pair (xi, Xj) according to the probability given the left hand side of formula (|3.2p can 
thus be done in O(N) operations rather than 0(N 2 ). All the required information 
can be held in binary tree structures (as described in [15]) whilst allowing an even 



9 This decrease is of the same order for the weighted particle system and for Marcus-Lushnikov's 
dynamics; it is worse for the coupled system. In this approach, at is approximated by the ratio 

\+lS\;N A-i<5A;iV „ A+i<5A;JV A-i<5A;iV 

(fi t — fj, t J/oA, where /j, t and /x t z are two coupled Markus-Lushmkov 

processes. So, the smaller <5A is, the more ^ + and /i^ (and /j, X+ ^ SX ' N and /i^ 2 SX ' N 

with it) look the same. This means that the 'real' number of particles in the difference fj, t +5 ' — 

/Ltj 2* A ' JV ; s a 'f unc ti on ' fsx(N) $J N of SX that decreases as SX goes to 0, a necessary condition 
for the ratio to be a good estimate of a±. 
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further reduction in the number of operations to choose each particle from 0(N) to 
O(logiV). Updating this information also requires O(logiV) operations. Further, the 
sums in the first fractions of the right hand side of p,2|) are automatically contained 
in the tree structure without further computation. 

Note that in the theoretical framework used in section the function <p{x) <p(y) 
can be used as a unique majorant kernel. We have yet chosen to present the above 
general procedure as we shall consider situations in which the above theory does not 
apply directly. 



3.3. Algorithm description. Recall 0^ is of the form ( jjXt, jyYt, jf^tj for 

a Markov process t = (X t , Y t ,Z t ) whose components are sums of Dirac masses and 
whose dynamics was described in section [27X1 What the algorithm really simulates is 
the discrete measure- valued process 04; a rescaling gives the time evolution of 0f . 
The algorithm is described in Algorithms Q] and [2] below. 

Note that there may be up to three different majorant kernels — for K, K + and 
K_. Therefore, we slice up the total majorant rates according to the event type 
a 6 {0, 1 + , 1~, 2+, 2~} to occur. We then have K a p such that J2a Ka/3 = K<* (from 
eq. l3.ip . where K a £ {K,K + ,K_}. This gives the corresponding rates p a p and p a . 

To order to describe numerical results it provides, we shall denote by L the number 
of simulations with the same initial conditions and by t run the computational time 
taken to run the algorithm (CPU time in seconds). 

4. Numerical Results. We have chosen to illustrate our approach in situa- 
tions where the theoretical results of section [2] do not apply, so as to show its ro- 
bustness. The main motivation of this article is to produce a stochastic estimate 
of the sensitivity a t whose variance is smaller than that given by existing methods. 
One step in this direction was done in [9], where cr t was approximated by the ratio 

(fit + 2SX ' N — /if 2(5A ' N )/SX, for two Marcus-Lushnikov processes with slightly dif- 
ferent parameters. The method there called for coupling them so as to reduce the 
variance of this estimator as much as can be done; this was done in the same spirit as 
the coupling used above. We shall refer to this algorithm as the CD algorithm (for 
central difference). The variance reduction obtained by this method is significant; 
we shall thus compare our results with those given by the CD algorithm. As our 
algorithm simulates a t directly, it will be called Exact; and depending on whether 
or not we use the coupling step we shall talk of the Exact Coupling or Exactlndep 
algorithm. 

The data presented deal with the additive kernel K(x, y) = X(x + y) and a kernel 
that is used in modelling soot formation in a free molecular regime [16} \T7\ [T8] , thus 
we shall call it the 'Soot Kernel': 



both are considered in the discrete setting where masses are integers. The reference 
value of A for the additive kernel will be 1 and for the soot kernel 2.1. We shall 
always take as initial condition for the Marcus-Lushnikov process N particles with 
mass equal to 1, and (Jq' N = (Tq' = 0. 
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Smoluchowski equation has an explicit analytic solution for an additive interaction 
(see the review by Aldous [19] for instance) we can compare our results with it; it will 
be convenient to write <7 t °° for at in this case. No analytic solution of Smoluchowski 
equation or its sensitivity equation is available for the soot kernel; we shall thus 
compare our estimators with what the ExactCoupling algorithm gives us for very 



Algorithm 1: The ExactCoupling algorithm - Part 1 



i Set t 

2 



0. while t < t en d do 



Generate a realisation of the holding time At with exponential law 

+ At. 

1~, 2 + , 2~} to occur with distribution 



of parameter J2 a p a , and set t <— t + At. 



Choose event type a € {0, 1 

Pa 

E„ Po • 

Choose process j3 with distribution fe^ . 

Given a and /3, choose a pair of particles using the index 

distribution 



Pa/3 



(3.3) 



where (xi,Xj) are the masses of particles sampled from the 
appropriate ensembles (fi N , a +,N or <r~ :N ) depending on a. 
Perform the coagulation step which depends on a: 
switch the value of a chosen do 

case a = 0; this part is the original Marcus- Lushnikov process. 
The chosen pair of particles is of the form (#,, Xj). 
With probability Jp make the jump 

AQ N = (S Xi+x . - ^ - S Xj ) ffi 0. 



case a = l + or 1 

The chosen pair of particles is of the form (xj, Xj). Set 

max{AT 2 + ,K 2 - } 



P 



and generate a realisation of a uniform 



random variable U in (0, 1). 
if U ^ p then 

if K 2 + > K 2 - then 

make the jump AQ N 
else 

make the jump AQ N 



0®5, 



Xi+Xj 



case a 



0®S Xi 
else Go to Step[TJ 

2+or 2"; Go to Algorithm H 



For each particle of a N that has just been involved in a coagulation 
or newly formed, do a cancellation operation if it can be done. 



16 STOP. 
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Algorithm 2: The ExactCoupling algorithm - Part 2 (Cases 
a = 2+,2" only) 

1 case a = 2 + or a = 2~ 

The chosen ordered pair of particles contains one particle of 
[i N and one particle of either order, 

if the pair is of the form (xi, •) where Xi is the mass of a 
particle from fi N then 

if the second particle belongs to <r + ' N then 
Choose a particle of a~' N according to the 
distribution 



10 



Efe[i,...,g] 9<*p{zi) 



(3.4) 



else 



Choose a particle of a + ' N according to the 
distribution 



Sfee[i,...,p] 9ap{.yk) 



(3.5) 



Set 

r+ ■= X] 9afi(yk) , 
fce[i,..., P ] 



9ai3{z e ) (3.6) 

ee[i,...,g] 



else 

Do the symmetrical operation, swapping g a( a with / a| g. 

The preceding steps produce a triple {xi,yk, zi) of particles 
from fx N ® a+' N ® a~ N . Set 

min{r + , r_} A" max{r_(_, r_} A' 



r+ + r_ K 



r+ + r- K 



Generate realisation of a uniform random variable U in 
(0,1). 

if < U ^ Pmin then 
make the jump 

A6 W = © ((y Xl+infc - y J © (<y X4+ „ - 5 Z< ). 
else if p min < U ^ p max then 
if r + > r_ then 

make the jump AQ N = (^j+3/ fc — <5y fc ) (5 Ki . 
else 

_ make the jump A@ N = ® S Xt ® (S Xi+Z( - S Ze ) . 



else 

_ Go to Step [D of Algorithm [TJ 



li Go to Step [U of Algorithm [TJ 
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01234 01234 
log(size) log(size) 



(a) ExactCoupling, t = 0.5 (b) CD (<5A = 0.05), t = 0.5 




log(size) log(size) 
(c) ExactCoupling, t = 3.0 (d) CD (<5A = 0.05), t = 3.0 

Figure 4.1. Sensitivity for additive kernel, X = 1.0, N = 10 3 , L = 1000. The confidence 
intervals for the larger particle sizes are omitted for pictorial clarity. 



high settings, say N = 3 x 10 6 and L = 10 3 simulations. Given any N, the I th run of 
the algorithm produces an estimator of a t which we shall denote by a t ' N . We shall 
set of° := 10 -3 E] a\' w . Figures |4~T1 and l4~2l show the empirical estimate of at 

i=l,...,10 3 

given after L runs, at different times. The line represents of°. For comparison, the 
results given by the CD algorithm for the same setting, with SX = 0.05, are plotted 
using stars. Also, Figure l4~3l shows what the solution to the original Smoluchowski 
equation looks like. 

To quantify the convergence of the empirical sensitivity 

-L-N 1 l-N 
1=1. .L 

to of° as N increases we have plotted in Figure l4~4l the quantity 

3 i>l 

where <f t L ' N (i) and represent the empirical and real sensitivities at particle 

mass i £ N respectively, and d vaT (N) represents the total variation distance between 
the empirical sensitivity and the sensitivity itself summed over some chosen time 
pointql^l {tj}. These results empirically confirm Theorem 12.21 (in this case where it 
does not apply), and quantify the speed of convergence as being of order jj. The 
analogue result for the CD algorithm is given in [9]. 



For Figure l4~4l the times points {tj} were chosen to be 0.125.J for j = 1, . . . , 56 
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O 



12 3 4 
log(size) 



12 3 4 
log(size) 



(a) ExactCoupling, t = 0.5 



(b) CD (SX = 0.05), t = 0.5 





12 3 
log(size) 



(c) ExactCoupling, t = 3.0 



(d) CD (SX = 0.05), t = 3.0 



Figure 4.2. Sensitivity for soot kernel, X = 2.1, N = 10 , L = 1000. The confidence intervals 
for the larger particle sizes are omitted for pictorial clarity. 



1 g 5 <r 
log(size) 



(5 1 2 3 4" 
log(size) 



(a) Additive kernel, t = 1.0 (b) Additive kernel, t = 3.0 



(5 1 2 5 4~ 
log(size) 



- o o o oo 

log(size) 



(c) Soot kernel, t = 1.0 



(d) Soot kernel, t = 3.0 



Figure 4.3. /xt as a function of log(particle size) 
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-1 slope 




4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 
log(N) 

(a) Additive kernel 




4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 
log(N) 

(b) Soot kernel 



Figure 4.4. Convergence in N of the ExactC oupling algorithm, N = 100 X 2 l for i = 0, 
NL = 2 x 10 8 . 



,5, 





Figure 4.5. log Varjv(i) as a function of N . The meaning of the symbols are as follows: Circles 
= ExactC oupling, Diamonds = Exactlndep, Triangles = C'D(5X = 0.10J, Crosses = C'D(5X = 0.05J, 
Pluses = CD(SX = 0.01). 



4.1. Variance. To analyse the variance of the random output of the algorithm 
we shall define the empirical variance at particle mass i £ N and time t as 



i=i 
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and shall take as a measure of the variance the quantity 

Varjv(t) —^Varivfrt). (4.1) 

Figure 14.51 represents its graph as a function of N using different algorithms. It 
shows that the ExactCoupling algorithm achieves a variance reduction by a factor 
10 3 compared to the CD algorithm. The plots also show that VarAr(t) is proportional 
to jj, a fact that should be related to a central limit theorem. 

4.2. Computational efficiency. Although section 14.11 indicates that the Ex- 
actCoupling algorithm produces very accurate estimators of the sensitivity, it comes 
at the price of a computational time greater than the one needed by the CD algo- 
rithm. This comes from the fact that the latter algorithm being essentially a Marcus- 
Lushnikov algorithm, it uses a generally decreasing amount of information, as the 
number of sensitivity particles decreases with time. On the other hand, the Ex- 
actCoupling algorithm has to deal with more and more sensitivity particles, whose 
number tends to grow exponentially. To see whether the gain of accuracy given by 
the ExactCoupling algorithm is worth the effort we propose two criteria. 

4.2.1. CPU time to reach a certain level of accuracy. Fix the observation 
time t (we choose large enough t so that the particle system has experienced many 
jumps, and therefore the variances are expected to be larger - see Figure l4~3l Given a 
certain level of accuracy, v, find for each algorithm the smallest A" for which Varjv(i) is 
smaller than v. See what computational time is needed to run the algorithm for this N 
(during an evolution time t for the particle system). Tables |4~T1 and [4721 show that the 
ExactCoupling algorithm remains mostly better than the CD algorithm. It also shows 
that it converges much quicker to the true sensitivity than the CD algorithm does. 
Note that for the soot kernel the CD algorithm with 5X — 0.1, 10 5 initial particles are 
not sufficient to reach the given level of accuracy; this setup already requires a CPU 
time equal to 1058.91 seconds. The comparison with the corresponding time for the 
ExactCoupling algorithm is greatly in favour of the latter. 

Table 4.1 
Additive kernel, v = 1.43 X 10~ 4 



t 


1.0 


1.0 


3.0 


3.0 




algorithm 

^run (s6Cs) 


ExactCoupling 

6500 

281.15 


CD (SX = 0.10) 

55000 

593.99 


ExactCoupling 

2100 

99.22 


CD («JA = 

16250 

213.34 


0.10) 






Table 4.2 
Soot kernel, v = 2.57 X 10~ 5 






t 


1.0 


1.0 


3.0 


3.0 




algorithm 

AT 

^run (s6Cs) 


ExactCoupling 

10000 

379.01 


CD (SX = 0.10) 

100000 
1058.91 

(v not reached) 


ExactCoupling 

6350 

382.15 


CD (SX = 

55000 

1104.24 


0.10) 
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4.2.2. Gain factor. Eibeck and Wagner introduced in [14] another quantity to 
compare the relative efficiency of two algorithms. Fix the observation time t. Given 
a setup (K(-, ■), N, L), denote by T EC (t) and T CB (t) the empirical mean CPU time 
needed by the ExactCoupling and CD algorithms to be run up to time t. Denote 
also by Var^ c (i) and Var^ s (i) the empirical variances given by formula l|4.ip when 
computed using ExactCoupling and the given algorithm 'alg' respectively. The gain 
factor of an algorithm over ExactCoupling, similar to that as introduced by Eibeck 
and Wagner, is defined here by the ratio 

T EC (t) Var% c (t) 
T**(t) Va4 s (t) 

It is related in some way to the analysis made in section [4.2.11 See section 5 of |14j . 
Figures 14.61 and 14.71 plot the reciprocal gain (its logarithm) as a function of time. 
Triangles, pluses and crosses represent data of the CD algorithm, for SX = 0.01, 0.05 
and 0.10 respectively, circles represent data of the Exactlndep algorithm, and the 
horizontal line at zero represents the threshold for ExactCoupling. 




01234567 1 01234567 1 01234567 

t t t 



(a) N = 10 3 (b) N = 10 4 (c) N = 10 5 



Figure 4.6. Additive kernel: log( Gain factor -1 ) as a function oft. The meanings of the 
symbols are as follows: Circles — Exactlndep, Crosses = CD(S\ = 0.10J, Pluses — CD(5X = 0.05 J, 
Triangles = C'D(8X = 0.01). The horizontal line is the threshold value 1.0 for ExactCoupling. 




Figure 4.7. Soot kernel: log(Gain factor -1 ) as a function oft. The meanings of the symbols 
are as follows: Circles = Exactlndep, Crosses = CD(8\ = 0.10J, Pluses = CD(8X = 0.05J, Triangles 
= CD(8X = 0.01J. The horizontal line is the threshold value 1.0 for ExactCoupling. 
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Figures l4~6l and l4~7l show good results. By and large, the CD algorithms appear 
to be considerably inferior to the ExactCoupling algorithm for the Soot kernel, and 
the Exactlndep algorithm in either performs slightly better than the CD (SX = 0.10). 
There appears to be little to moderate difference in behaviour over different values of 
N. 

The picture is different for the Additive kernel. For N = 1000, we find that the 
CD (<5A = 0.1) is better than the ExactCoupling, at least for very small or large times. 
This disadvantage gradually disappears over larger N — this is due to the increased 
probability of cancellations for larger N which reduces the number of particles in 
the ensembles and therefore the CPU times. Other than this, the ExactCoupling 
algorithm maintains a substantial lead over the other algorithms. 

5. Conclusions. A stochastic particle system approximation to the parametric 
sensitivity in Smoluchowski's coagulation equation was introduced. Rather than tak- 
ing a finite difference approach to calculating sensitivities, we considered the direct 
parametric derivative of (|1.2p . and developed a Monte-Carlo algorithm which would 
approximate its solution. The particle system approximation was proved to converge 
weakly to the solution of the sensitivity equation (|1.2p . as the number of particles 
increases indefinitely. 

The first algorithm developed (Exactlndep) allows for an exponential increase in 
the number of sensitivity particles. We sought to reduce this increase using several 
tricks: Cancellation removes 'unnecessary' sensitivity particles which are needed to 
describe it, whilst coupling prevents their creation. These make a significant reduction 
to the number of particles in the ensemble. Furthermore, the resampling method puts 
a cap on the total number of sensitivity particles, thus stopping their exponential 
escalation. This gives us the ExactCoupling algorithm. 

In the Numerical Results section, it was empirically confirmed that the order of 
convergence is 0(1/ N) where N is the number of initial particles. We then compared 
the Exact algorithms with those found in [9j, named here CD algorithms. It was shown 
that the variance of the sensitivity estimators were orders of magnitude smaller for 
the ExactCoupling algorithm than for the CD algorithms. However this came at the 
price of longer CPU run times. Two measures of efficiency, taking both the variance 
and the CPU time into account, were then considered. The ExactCoupling algorithm 
happens to require much smaller time to to reach a fixed level of error than any CD 
algorithm, and the gain factor, as defined in [14], also happens to be in favour of the 
ExactCoupling algorithm, most of the time. This definitely gives a clear advantage 
of our approach over finite difference methods. 

However, both methods have some inherent drawback: unlike the adjoint method 
[8] , they are unidimensional in nature and compute sensitivity only for a fixed value of 
the parameter. It would be useful to construct a particle system approximation which 
do not have these weaknesses. Also, although the convergence theorem established 
in section [2] in a general framework is quite encouraging, it is not clear whether the 
algorithm will be as efficient as above if particles's masses can take any positive value. 
We leave the investigation of these questions for future work. 
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